install.packages('sf')
install.packages('ggplot2')
install.packages('tidyverse') 
install.packages('gridExtra') 
install.packages('patchwork')
install.packages('readxl')
install.packages('haven')

rm(list = ls())


library(sf)
library(ggplot2)
library(tidyverse)
library(gridExtra)
library(patchwork)
library(readxl)
library(haven)


#setwd("~/Dropbox/05 Bismarck/05 Map/")
setwd("C:/Users/Felix Kersting/Dropbox/05 Bismarck/05 Map")

###############################################################################################################
###############################################################################################################

# read matching file
match <- read_excel("mapid_eid.xlsx") %>% filter(LAND == 1000)

# read shapefile
shape <- st_read("German_Empire_1871_v.1.0.shp") %>% select(ID) 

# merge, need to put shapefile first, so that geometry is kept
sm <- merge(shape, match, by='ID')

# aggregate ny e_id, only keeps e_id and geometry
# 229 obs
sm_agg <- sm %>% 
          group_by(bismarck_id) %>% 
          summarise() 

# plot to see
sm_agg %>% ggplot() + geom_sf()

# write as shapefile to use for maps (after merging with data by e_id)
st_write(sm_agg, "bismarckmap.shp", append=FALSE)

# centroids, transform to wsg84 (EPSG 4326)
points_xy <- sm_agg %>% st_centroid() %>% st_transform(4326) %>% mutate(X = st_coordinates(.)[, 1], 
                                                                   Y = st_coordinates(.)[, 2])

# plot to see
points_xy %>% ggplot() + geom_sf()

# write as dta-file to use with acreg
points_xy %>% st_drop_geometry() %>% write_dta(., "bismarckmap.dta")


mapdata <- read_excel("map_data.xlsx") 

#shape <- st_read("bismarckmap.shp") 

gdf <- merge(sm_agg, mapdata, by='bismarck_id')

#Map 1
gdf <- gdf %>% mutate(quantiles = as.factor(ntile(s_add_treat, 10)))
p <- gdf %>% ggplot() + geom_sf(aes(fill = quantiles), size = 0.1)  + 
scale_fill_brewer(palette = "RdYlBu", name = "Newly-insured workers", direction = -1, guide = guide_legend(nrow=2)) + 
theme_void() + theme(legend.position = "bottom")
p

pdf("fig_reform.pdf")
plot(p)
dev.off()

#Map 2
affected_breaks <- c(0, 0.5, 1.1)
affected_labels <- c("Not affected", "Affected")

gdf <- gdf %>% mutate(d_socialist = cut(d_verboten, breaks = affected_breaks, labels = affected_labels, right = FALSE))

p <- gdf %>%ggplot() +
  geom_sf(data = gdf, aes(fill = d_socialist)) +
  scale_fill_manual(values = c("white", "red"), name = "Anti-socialist law") + theme_void() + theme(legend.position = "bottom")
p

pdf("fig_anti.pdf")
plot(p)
dev.off()

#Map 3
affected_breaks <- c(0, 0.5, 1.1)
affected_labels <- c("No readers", "Readers")

gdf <- gdf %>% mutate(d_radical = cut(d_sozialdemokrat, breaks = affected_breaks, labels = affected_labels, right = FALSE))

p <- gdf %>%ggplot() +
  geom_sf(data = gdf, aes(fill = d_radical)) +
  scale_fill_manual(values = c("white", "red"), name = "Readers forbidden newspaper") + theme_void() + theme(legend.position = "bottom")
p

pdf("fig_radical.pdf")
plot(p)
dev.off()
